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Bursts of action potentials within neurons and throughout networks are believed to serve 
roles in how neurons handle and store information, both in vivo and in vitro. Accurate 
detection of burst occurrences and durations are therefore crucial for many studies. A 
number of algorithms have been proposed to do so, but a standard method has not been 
adopted. This is due, in part, to many algorithms requiring the adjustment of multiple 
ad-hoc parameters and further post-lioc criteria in order to produce satisfactory results. 
Here, we broadly catalog existing approaches and present a new approach requiring the 
selection of only a single parameter: the number of spikes N comprising the smallest 
burst to consider. A burst was identified if N spikes occurred in less than T ms, where 
the threshold T was automatically determined from observing a probability distribution of 
inter-spike-intervals. Performance was compared vs. different classes of detectors on data 
gathered from in vitro neuronal networks grown over microelectrode arrays. Our approach 
offered a number of useful features including: a simple implementation, no need for ad-lioc 
or post-hoc criteria, and precise assignment of burst boundary time points. Unlike existing 
approaches, detection was not biased toward larger bursts, allowing identification and 
analysis of a greater range of neuronal and network dynamics. 
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INTRODUCTION 

Bursts of action potentials within neurons and throughout net- 
works are believed to serve roles in how neurons handle and 
store information, both in vivo and in vitro. Accurate detection of 
burst occurrences, durations, and boundaries are therefore cru- 
cial for many studies. In vivo studies have linked bursting, also 
referred to as "up states" in some cases, to oscillations travel- 
ing through brain regions, encoding of sensory input, assisting 
transmission of neural information, and markers of disease states 
such as epilepsy (Lisman, 1997; Izhikevich et al., 2003; Staley and 
Dudek, 2006). The prominent feature of primary central nervous 
system cells, cultured in vitro, is global network bursts, which are 
interspersed by tonic activity of varying degrees within a portion 
of neurons. Burst statistics have been used, for example, to study 
how information could become encoded after applying electrical 
(Madhavan et al., 2007) or chemical (Eytan et al., 2004; Selinger 
et al, 2004) stimulation to induce network plasticity. They have 
also been used to judge in vitro developmental stage, whereby 
networks transition from sparse uncorrelated spiking until dis- 
playing synchronized and aperiodic network bursts of various 
magnitudes in a "mature" network after 3-4 weeks (Van Pelt et al, 
2004; Wagenaar et al., 2006). Both in vivo and in vitro, finite and 
repeating sets of activity, variously termed motifs or songs or 
assemblies, have been identified and proposed to be substrates to 
store memory traces (Baruchi and Ben-Jacob, 2004; Ikegaya et al. 



2004; Segev et al., 2004; Harris, 2005; Eytan and Marom, 2006; 
Rolston et al, 2007; Kumar et al, 2010). 

While a standard method to identify and detect network bursts 
has not been adopted, a number of algorithms were proposed 
for both in vivo (Legendy and Salcman, 1985; Cocatre-ZUgien 
and Delcomyn, 1992; Kaneoke and Vitek, 1996; Elias et al, 2007; 
Gourevitch and Eggermont, 2007; Ji and Wilson, 2007; Ko et al., 
2012) and in vitro spike trains (Mukai et al., 2003; Xia et al., 
2003; Segev et al, 2004; TurnbuU et al, 2005; Wagenaar et al, 
2005; Selinger et al, 2007; Pasquale et al, 2010; Tokdar et al, 
2010; Pimashkin et al, 2011; Kapucu et al, 2012; Weihberger 
et al., 2013). Our working definition of a burst will be a period 
of high-frequency occurrences of multiple action potentials inter- 
spersed by periods of lower frequency tonic activity. We consider 
a network hurst as synchronous spikes spatially distributed across 
multiple recording channels, including cases when few spikes 
occurred per channel. Burst detectors aim to build a border 
separating the higher and lower activity regimes. A difficulty 
arises because clear borders are not necessarily apparent and 
vary between, and even within, preparations. Many detection 
algorithms therefore require adjusting multiple ad-hoc param- 
eters and farther post-hoc criteria in order to produce visually 
acceptable results. 

The various burst detection algorithms can be broadly cat- 
egorized into two classes: (1) those setting a rate-threshold to 



Frontiers in Computational Neuroscience 



www.frontlersin.org 



January 2014 | Volume 7 | Article 193 | 1 



Bakkumet al. 



Parameters for burst detection 



detect bursts whenever the activity rate exceeds a specific value; 
(2) those setting an inter-spike-interval or ISI-threshold to detect 
bursts whenever the ISI between consecutive spikes is less than 
a specific value. Rate-threshold detectors simply bin together 
the spike times from all recording channels within a specified 
time window in order to create a firing rate histogram. In its 
most basic implementation, two parameters need to be set: the 
time window and the activity rate threshold. A bursting regime 
is then identified whenever the number of spikes exceeds the 
threshold (Mukai et al, 2003; Xia et al, 2003; Ji and Wilson, 
2007; Pimashkin et al., 2011) or, similarly, whenever the num- 
ber of active electrodes exceeds the threshold (Segev et al., 2004). 
This method works better for detecting multi-channel network 
bursts. This is due to the fact that network-wide spike trains 
provide higher signal-to-noise ratios than single-channel spike 
trains: spike counts from merged single-channel spike trains sum- 
mate to multiplicatively larger values in time windows including 
bursts than those during non-bursting periods. The parameters 
are usually chosen empirically, in part, because the distribution 
of histogram peaks is often continuous. However, for a given time 
window, a rate threshold can be automatically set as the spike (or 
electrode) count value that separates the peaks in a bimodal prob- 
ability distribution of firing rates ( Ji and Wilson, 2007), also called 
a "discharge density" (Kaneoke and Vitek, 1996). Burst bound- 
ary time points are approximated at detection threshold crossings 
or via additional threshold settings. Histogram peaks and thresh- 
old crossings have been more precisely found by convolving spike 
times with Gaussian kernels, a decay function, or other smooth- 
ing methods (Xia et al., 2003; Segev et al., 2004; Ji and Wilson, 
2007). The tendency to choose "safe" thresholds favors the detec- 
tion of larger bursts, which underestimates burst number and 
duration. 

ISI-threshold detectors consider that periods of low and high 
ISIs correspond to spikes occurring within and outside of bursts, 
respectively, and analyzing peaks in the probability distribution of 
the ISIs can identify appropriate thresholds (Cocatre-Zilgien and 
Delcomyn, 1992). The ability to assign boundary time points to 
specific spikes offers an advantage over rate-threshold detectors. 
At its most basic implementation, the ISI threshold is the only 
parameter required and can be automatically selected (Pasquale 
et al, 2010). A number of simple and complex methods to guide 
ISI threshold selection have been proposed. These are commonly 
based on finding valleys in the distributions of plain ISIs or 
the logarithm of ISIs or from the discharge density (Kaneoke 
and Vitek, 1996; Wagenaar et al, 2005; Selinger et al, 2007; 
Pasquale et al., 2010; Kapucu et al, 2012). The latter, actually, 
transforms the spike count value separating bimodal probabil- 
ity distributions of firing rates, as opposed to that of ISIs, to 
an ISI threshold by multiplying its inverse by the time window 
used to calculate the firing rates. Often however, a number of 
post-hoc criteria are introduced to better fit the data, includ- 
ing reintroducing rate-based metrics such as minimum spike 
counts or number of channels activated. A separate important 
branch of ISI threshold detectors statistically compares recorded 
ISIs to what would be expected assuming spike activity behaved 
following a model distribution. Burst regimes are then identi- 
fied whenever activity exceeds expectations, termed a "surprise." 



The famous Poisson Surprise method was introduced many 
decades ago (Legendy and Salcman, 1985) and amended variously 
into non-parametric Rank Surprise (Gourevitch and Eggermont, 
2007), Robust Gaussian Surprise (Ko et al., 2012), and Pause 
Surprise (Elias et al., 2007) methods. Unlike rate-threshold detec- 
tors, ISI-threshold detectors typically operate on single-channel 
spike trains. Single-channel burst events detected in the first-stage 
are then combined together in order to identify network bursts, a 
process requiring additional parameters (Wagenaar et al, 2005; 
Pasquale et al., 2010). A main reason for the first-stage detection 
is because cumulative ISI distributions from multiple channels 
tend to average out peaks that would be apparent in probabil- 
ity distribution plots of single-channel spike trains. This obscures 
the choice of ISI threshold. Any method using a first-stage iden- 
tification of single-channel bursts will also be biased toward the 
detection of large network bursts: network bursts composed of 
synchronized spikes, but with only one or a few spikes per chan- 
nel, will be missed during the first-stage. Smaller "spatial" bursts 
of action potentials are reminiscent of the concepts of neuronal 
assemblies and synfire chains (Abeles et al., 1994; Kumar et al, 
2010), and their detection may therefore be useful. 

Our goal was to develop a simple yet robust network burst 
detector that does not depend on ad-hoc or post-hoc detection 
criteria and is able to detect smaller network bursts. Therefore, 
we built an ISIm -threshold detector, where ISIm is the inter-spike- 
interval between every N*^ spike instead of every consecutive 
spike. The key simple consideration is the fact that the time points 
of a set of N consecutive spikes will give a better representation 
of a network's firing rate status than 2 consecutive spikes (i.e., 
ISI). Only one detection parameter is set: the number of spikes N 
that compose the smallest network burst to consider. N = 10 was 
used throughout this paper and corresponds to approximately 
0.1 spikes per recording channel for our setup (126 channels). 
However, ISIn -threshold detection can be performed for a wide 
range of N for both single-channel and cumulative network spike 
trains. Peaks in the probability distribution of ISIm represent- 
ing bursting and non-bursting regimes become more and more 
apparent for increasing N. The border between burst regimes is 
again automatically assigned at the valley between the (logarith- 
mic) ISIn peaks. Compared to an ISI threshold, an ISIjv threshold 
will be easier to identify, produce fewer false-positive detections, 
maintain precise assignment of burst boundary time points, and 
allow simplification of the detection algorithm. Post-hoc ver- 
ification conditions are not required, and, like rate-threshold 
detectors, detection can be performed directly on the cumula- 
tive network spike train. ISIn -threshold detection is notated as: 
if r;_|-(j^_i)-r, <ISIn -threshold, then spikes S,- to S,-|-(n-i) are in 
the same burst, where T, is the time of spike S;; Matlab code is 
provided as Supplementary Material. 

We applied the ISI^ -threshold detector to data gathered from 
cultures of rat primary neurons and glia grown over complemen- 
tary metal-oxide-semiconductor (CMOS) -based microelectrode 
arrays (MEA), although the specific choice of recording device 
was not critical. From our judgment (see Discussion), the detec- 
tor performed well and avoided biasing identification toward 
larger bursts. In many cultures, a clear distinction between large 
bursts across the majority of channels and smaller bursts across 
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a subset of channels was apparent from observing distributions 
of the number of spikes or number of contributing channels in a 
burst. On the other hand, a continuum of smaller network bursts 
existed without a clear cutoff between bursting and non-bursting 
regimes. Such bursts would not be reliably identified with existing 
detectors, and much diversity in network information processing 
may be overlooked if smaller bursts are not considered. While the 
ISIj/ -threshold detector performed well by setting only a single 
parameter, additional criteria can improve performance. In par- 
ticular, large amounts of intermittent tonic spiking will degrade 
the performance of most detectors, and a method to identify 
and exclude such channels in an optional preprocessing step is 
presented. 

METHODS 
CELL CULTURING 

Techniques have been developed to maintain neural cultures 
and conduct experiments for many months (Hales et al., 2010; 
Bakkum et al., 2013). Briefly, E18 Wistar rat cortices were 
dissociated using trypsin and mechanical trituration. 20k-40k 
neurons and glia were seeded over an area of ~12mm^ on 
top of the CMOS chip. Layers of poly(ethyleneimine) followed 
by laminin were used to adhere cells. Plating media consisted 



of Neurobasal-B27 supplemented with 10% horse serum and 
0.5 mM GlutaMAX during the first 24 h. Growth media con- 
sisted of DMEM supplemented with 10% horse serum, 0.5 mM 
GlutaMAX, and 1 mM sodium pyruvate. Cultures matured for 
about 1 month prior to experimentation, and experiments were 
conducted inside an incubator to control of environmental con- 
ditions (36°C and 5% CO2). Burst detection was performed on 
9 cultures from 4 platings. Typical network activity is depicted in 
Figure 1. 

CMOS-BASED MEA AND RECORDING OF NETWORK ACTIVITY 

Cortical networks were grown for many weeks over 11,011- 
electrode CMOS-based MEAs (Prey et al, 2010; Livi et al, 
2010), which provide enough spatial and temporal resolution 
to detect action potentials from any neuron lying on the array: 
1.8 X 2.0 mm^ area containing 8.2 x 5.8 |xm^ electrodes with 
17.8 |i,m pitch (3150 electrodes per mm^), sampled at 20 kHz. 
Subsets of 126 electrodes can be read-out (and/or stimulated) 
at one time, and electrode selection can be re-configured within 
a few ms. Custom software on a personal computer, includ- 
ing modified Meabench code, a field-programmable gate array 
(PPGA), and a microcontroller embedded in a custom circuit 
board were used to acquire data (Wagenaar et al, 2005; MuUer 
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FIGURE 1 I (A) Network bursts (indicated in red) of varying durations (channels with gray dots). (B) Locations of the recording electrodes 

of coincidental APs (raster dots) recorded across multiple channels. selected in the CMOS-based MEA (see Methods). Gray triangles 

The majority of neurons mainly fired together within network bursts correspond to the channels with gray dots in (A). Scale bar: 

(columns of dots), but some also fired "tonically" outside of bursts 200 |im. 
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et al., 2013). To identify the locations of neurons growing over 
the array, a sequence of about one hundred recording config- 
urations were scanned across the whole array while recording 
spontaneous activity. To sample network activity, 126 recording 
electrodes were arbitrarily selected. Because an action potential 
from a single soma can be detected on multiple nearby elec- 
trodes, caution is required if configured electrodes are close to 
each other: activity from a single soma could be falsely detected 
as a small network burst. This case was avoided by maintaining 
a minimum inter-recording-electrode distance exceeding the spa- 
tial spread of somatic signals (about 40|xm). Alternatively, this 
case could be avoided by transforming channel spike trains into 
neuronal spike trains via spike sorting techniques (Franke et al., 
2012a,b; Jackel et al, 2012). Matlab R2012a was used for data 
analysis. 



PROPOSED NETWORK BURST DETECTION ALGORITHM 

The proposed ISIn -threshold method to detect bursts is depicted 
in Figure 2, and Matlab code is provided as Supplementary 
Material. To summarize, bursting (low ISI^r) and non-bursting 
regimes (high ISIn) typically formed two peaks in a histogram 
of the base-10 logarithm of ISIn. ISIn is the inter-spike-interval 
between every Nth spike in the network. The time points from 
all spikes on all channels were combined into a single train and 
input to the detector; spatial information about electrode loca- 
tions was not incorporated. Negative signals exceeding 5 standard 
deviations of the noise were considered to be spikes arising from 
somatic action potentials. The border between regimes, or ISI^; 
threshold, was chosen as the valley between the peaks (or the first 
minima for multiple peaks). Setting histogram bin widths to be 
equally spaced on a logarithmic scale increases the height of the 




FIGURE 2 I ISI/v thresholds and burst detection. (A) The probability of 
elapsed times between consecutive spikes (black; ISI) and every A/th spike 
(gray) up to every 10th spike (red) are plotted. Elevated firing during network 
bursting corresponds to lower ISI, and the red arrow indicates the threshold 
for burst detection used in Figure 1. (B) The elapsed time between 
consecutive spikes is plotted vs. the elapsed time between every 10 spikes. 
Histograms correspond to the black and red probability distributions in (A), 
and the red and black arrows correspond to the ISI and ISl^ thresholds in (A). 
For (A) and (B), ISIs were jittered by a random value between zero and one 



sample (50 (is) in order to better visualize the contribution from low ISIs. 
These would otherwise be plotted on top of each other in discrete lines 
corresponding to multiples of the sampling rate. The inset pie chart shows 
the percentage of spikes in each quadrant. Symbols match the spike markers 
in (C). (C) Detector performance for a segment of network activity. Black 
pluses and blue squares indicate spikes that would be classified in bursts (red 
bars) using an ISI/v^io threshold. Black pluses and cyan circles indicate spikes 
that would be classified in bursts according to an ISIn^2 threshold. Green 
triangles indicate spikes outside of bursts for either case. 
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non-bursting (high ISIjv) peak, providing better discrimination 
(Selinger et al., 2007; Pasquale et al., 2010). The user must choose 
a single parameter: the number of spikes N that compose the 
smallest network burst to consider. Detection is simply: If N con- 
secutive spikes occur within a time period equal to or less than the 
ISIn threshold, the spikes are assigned to a burst. The burst ends 
when this condition is no longer met. Since all spikes are assigned 
to be within or outside of a burst, specific time points for the first 
and last spikes in a burst are assigned. For this paper, we chose 
iV = 10 spikes, and identified ISIjv thresholds from 1 h record- 
ings. Channels with high levels of tonic spiking were excluded 
prior to analysis as depicted in Figure 5. This corresponded to 
between 2 and 10 percent of the recording channels (indicated 
by gray markers in raster plots). 

EXISTING BURST DETECTION ALGORITHMS USED FOR COMPARISON 

Rate-threshold detector 

A rate-threshold detector algorithm was adapted from the work of 
Ji and Wilson (2007). Specifically, a network burst was detected 
if a firing rate histogram with 50 ms or 5 ms time windows 
exceeded N spikes. Rate-thresholds, N, were automatically set 
at the valley between peaks in the spike count (or electrode 
count) probability distributions. A 50 ms time window corre- 
sponds to the choice by Mukai et al. (2003), who used a sim- 
ilar preparation to ours. A shorter time window (e.g., 5 ms) 
may allow more precise identification of burst start and end 
times. 

ISI-threshold detector 

An ISI-threshold detector algorithm was adapted from pseudo 
code provided in a recent and thorough publication by Pasquale 
et al. (2010), who also used a similar preparation. In the first-stage 
single-channel burst detection, ISI thresholds were automatically 
set at a valley between peaks in the (logarithmic) ISI probability 
distribution and capped at maximum value of 100 ms. A single- 
channel burst was detected if the ISIs of 5 consecutive spikes were 
each less than the ISI threshold. Network bursts, or a "burst of 
bursts," were detected using the same algorithm but on burst 
events instead of spikes. A post-hoc criterion of a minimum num- 
ber of activated channels per network burst was specified as 20% 
of the recording channels. This corresponds to 25 channels in our 
array, but a value of 10 channels was used to improve consistency 
with the ISIjv=io-threshold detector. 

Rank surprise detector 

The non-parametric Rank Surprise algorithm was acquired from 
Matlab code provided with the original publication by Gourevitch 
and Eggermont (2007). As data in the original paper were 
recorded in the cortex of an anesthetized cat in vivo, two default 
parameters were changed to improve consistency with the other 
detectors: The minimum number of spikes and maximum ISI 
within a burst were changed to 5 spikes and 100 ms, respectively. 
However, the default values of 3 spikes and the 75th percentile 
of ISI produced similar results. A method to identify network 
bursts was not provided. Therefore, we identified a network burst 
simply whenever at least 10 single channel bursts overlapped 
temporally. 



RESULTS 

Bursts were identified in 1-h recordings of spontaneous activity 
from primary cortical cultures using the ISIn -threshold detec- 
tor with N = 10 (see Methods). Figure 1 shows typical neuronal 
network activity. Bursts showed a variety of sizes and inter- 
burst-intervals. Tonic activity commonly occurred in a subset 
of channels to varying degrees. Some channels showed elevated 
tonic activity throughout the duration of a recording, while oth- 
ers showed intermittent tonic activity. Recording channels were 
ordered by the overall firing rate in order to better observe the 
occurrence of bursts, indicated by the columns of spikes in the 
raster plot. 




12 3 

10 10 10 

Network burst size [spikes] 

FIGURE 3 I Distribution of burst sizes. Increasing burst size correlated to 
increasing number of contributing channels (A) and increasing burst width 
(B). Larger bursts covering the majority of channels and smaller bursts 
covering a subset of channels are clearly distinguishable by observing 
valleys in each of the three histograms. 
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The method to find the ISIn=io threshold used to detect 
bursts in Figure 1 is depicted in Figure 2. The deflection in 
the probability curve suggests the existence of more than one 
type of activity state, such as bursting (low ISl^) and non- 
bursting (high ISIn; contributions from tonic activity) regimes. 
Increasing N increased the amount of deflection in the curve 
(Figures 2A, 4). This will add precision when determining a 
threshold and improve performance at the expense of increas- 
ing the criterion for the minimum number of spikes per burst. 
The ISIjv=io threshold was 50 ms for the represented data and 
150 ± 100 ms (mean ± SD) across 9 cultures. For compari- 
son, an ISIn=2 threshold was approximately 20 ms (Figure 2A, 
large black arrow) and 35 ± 45 ms (mean ± SD) across 9 cul- 
tures. The probability distributions and thresholds will depend 
on the overall number of neurons and their firing profiles, which 
may vary in time. The additional spikes that would be classi- 
fied as composing a burst for ISI]v=2 but not for ISI]v=io are 
presented as cyan circles in Figures 2B,C. While some of these 
spikes appear to be false positives, the close alignment of others 
may represent bursts, or perhaps assemblies, with fewer than 10 
spikes. 

In 6 out of 9 cultures, a clear distinction between large bursts 
covering the majority of channels and smaller bursts covering 
a subset of channels was apparent from observing distributions 
of burst size (number of spikes or elapsed time) or the num- 
ber of contributing channels (Figures). Small bursts averaged 
less than 1 spike per recording channel, indicating they would 
be difficult to identify using existing detectors. A continuum of 
smaller bursts existed without a clear cutoff between bursting and 
non-bursting regimes. This suggests using caution when mak- 
ing assumptions about burst dimensions or "surprise" criteria for 
burst detection. Considering the small bursts represented 80 per- 
cent of all bursts and 10 percent of the spiking for this culture, 
much diversity in network information processing may be over- 
looked if only large bursts are detected and analyzed. Small bursts 



in the other 5 cultures represented 76, 45, 26, 25, and 12% of 
all bursts for each culture. From an information theory perspec- 
tive, the extreme cases of no neuron or all neurons firing (i.e., the 
largest bursts) provide no information content, or only 1 bit if 
combined as on and off states. Smaller bursts then theoretically 
hold more information content. A previous experiment demon- 
strates the utility of small bursts, where electrically evoked bursts 
were used to instruct motor output in an embodied cultured 
network (Bakkum et al, 2008). Successful goal-directed behav- 
ior, based on plasticity in the spatio-temporal burst structure, 
was possible only when using a stimulation electrode that evoked 
small bursts. 

While the ISIjv-threshold detector was robust in the pres- 
ence of tonic activity, which in fact helped to form the sec- 
ond peak and valley in the logarithmic ISI^ histograms that 
were used to determine the ISIm threshold (Figure 4), strong 
tonic firing on too many channels may compromise detection 
performance. For example, sustained tonic spiking could cause 
neighboring bursts to be falsely identified as a single burst. 
Also, tonic spiking on multiple channels could produce false- 
positive bursts. In any case, channels exhibiting high tonic activity 
are readily apparent from a visual inspection of raster plots 
(Figure 1). As an optional pre-processing step, these channels 
can be excluded. To make exclusion less subjective and more 
automated, a method to quantify the amount of tonic activ- 
ity on a given recording channel is depicted in Figure 5. Here, 
a "tonic firing rate" was calculated by enforcing an arbitrarily 
chosen refractory period of 250 ms, whereby any spike follow- 
ing within the refractory period was excluded. This filtered out 
most spikes occurring during bursts while preserving much of the 
tonic spiking occurring between bursts. Channels with high tonic 
firing rates corresponded well to visual observations (Figure 5). 
A suitable choice of refractory period is longer than the aver- 
age ISI of a neuron and shorter than the average inter-burst- 
interval. 



Excluding tonic chianneis 



Including tonic channels 




FIGURE 4 I Valleys in the bimodal ISI/v probability distributions deepen 
with increasing N (light gray to dark gray) or with the inclusion of 
channels exhibiting tonic spiking (right panel). The distributions (black, 
red, light gray to dark gray) correspond to W equal to 2, 10, 50, 100, 250, 500, 



threshold, approaches the maximum burst size (~4000 spikes), which 
represents the trivial case of no spikes being within a burst. 8 out of 102 
channels were identified as tonic according to the method presented in 
Figure 5. The inclusion of tonic channels can improve valley identification with 



1000, 2000, and 4000. The first peak disappears as W, the minimum burst size a tradeoff of increased risk of identifying consecutive bursts as a single burst. 
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5 10 15 

Total Firing Rate [Hz] 



FIGURE 5 I Quantifying levels of tonic activity. (A) Example network 
activity exiiibiting bursting (red bars) and cliannels with high levels of tonic 
spiking (gray triangles). Spikes on all other active channels are indicated by 
black pluses. (B) The same network activity after introducing an artificial 



refractory period of 250 ms. Spikes during bursting became preferentially 
removed. (C) For each channel, the tonic (from B) vs. total (from A) firing rate 
is plotted. Channels with elevated tonic firing (gray triangles) match visual 
observations in (A) and (B). 



DISCUSSION 

COMPARISON TO EXISTING DETECTORS 

Figures 6, 7 provide a comparison of the ISI^ -threshold detec- 
tor with detectors representative of the rate -threshold, ISI- 
threshold, and surprise methods by Ji and Wilson, Pasquale 
et al., and Gourevitch and Eggermont, respectively (Mukai 
et al, 2003; Gourevitch and Eggermont, 2007; Ji and Wilson, 
2007; Pasquale et al, 2010). WhUe the rate-threshold and ISI^- 
threshold methods operated directly on the cumulative network 
spike train, the ISI-threshold and the Rank Surprise meth- 
ods required a first-stage single-channel burst detection step 
prior to identification of a network burst (Figures 6B,C). The 



ISIn=io threshold was equal to 140 ms in Figure 6 and 50, 250, 
and 350 ms in Figures 7A-C, respectively. The rate threshold 
for a 50 ms time window was equal to 8 spikes in Figure 6 
and 160, 80, and 10 spikes in Figures 7A-C, respectively. Rate- 
threshold detector performance was similar for a 5 ms time 
window (see Figures). Refer to the Methods for implemen- 
tation details. To compare performance in different culture 
preparations, Figure 7 contains data from 3 cultures whose 
small bursts were, according to the ISIjv=io-threshold detec- 
tor, 80, 12, or 0% of the total number of bursts per each 
culture. Data from the culture in Figure 7 A is also presented 
in Figures 1-3. 
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FIGURE 6 I Performance of the ISI/v=io-threshold detector compared 
to example rate-threshold, ISI-threshold, and Rank Surprise 
detectors. (A) Burst times and durations for each detector (color bars 
and legend) applied to the network activity presented in the raster plot 
(each marker is one spike). The ISI/v=io-threshold and rate-threshold 



detectors were performed directly on the network spike train. The 
ISI-threshold and Rank Surprise detectors were performed on 
single-channel spike trains in a first stage. The respective single-channel 
burst events in (B) and (C) were then combined in a second stage to 
identify network bursts as described in the Methods. 



All methods detected the largest network bursts, and the 
ISIn -threshold detector also detected smaller network bursts 
(Figures 7 A-C; left column). In Figure 7A, 20% of the bursts 
detected by the ISIjv=io -threshold method were considered to 
be large bursts (c.f Figures). This is consistent with the ISI- 
threshold, rate-threshold, and Rank Surprise methods detecting 
18, 16, and 11% of the total bursts identified by the ISIn=io- 
threshold method, respectively. For Figure 7B, the values were: 
88% large bursts; 90, 90, and 76% of the total detected. For 
Figure 7C, the values were: 100% large bursts; 100, 100, and 
93% of the total detected. The ISI-threshold detector strug- 
gled with tonic activity, whereby some bursts had artificially 
extended durations that sometimes merged neighbored bursts 
(Figure 6). The authors of the ISI-threshold detector acknowl- 
edged this situation (Pasquale et al., 2010), and they and others 
provided a post-hoc solution described below. Merged bursts 
are also possible with the ISIn -threshold detector, especially 
when the ISIn threshold is relatively large. The rate-threshold 



detector performed well. However, burst durations were con- 
sistently underestimated (Figures 7A-C; right column), and the 
smallest bursts in a more challenging dataset were not detected 
(Figure 8, Supplemental Figure 1). Furthermore, the inability to 
accurately assign burst boundaries can be a drawback for some 
applications. Interestingly, rate-threshold detection with a shorter 
5 ms time window underestimated burst duration to a greater 
extent (Figure 8). A rate-threshold detector based on the number 
of electrodes or the number of (spike-sorted) neurons displaying 
activity, as opposed to the number of spikes recorded, produced 
similar results (Figure 8, Supplemental Figure 1). Rate-threshold 
detector performance will improve with increasing numbers of 
recording channels: time windows that contain bursts will sum- 
mate to multiplicatively larger values than windows containing 
lower frequency tonic activity. 

The Rank Surprise method detected the fewest bursts and most 
underestimated burst duration (Figure 7). A commonly touted 
advantage of surprise methods over other methods is their ability 
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FIGURE 7 I An ISI/v=io-threshold detector identified small-sized bursts, 
and measured burst durations were shorter for a firing-rate-threshold or 
a rank surprise detector (A-C) Each method was applied to 1-h long 
recordings from 3 cultures that each had different amounts of small bursts. 
The percentage of small bursts out of the total number of bursts, according to 
the ISI)v=io-threshold detector, were 80, 12, and 0% for (A-C), respectively. 
Each detector identified the largest bursts, and the ISIw^io-threshold 
detector identified smaller bursts (left column). Since different detectors 
assign different durations to the same burst, one detector (ISIw^io) was 



chosen as a reference (x-axis). In this manner, the same burst will be plotted 
at the same x-location, and its occurrence can be compared across detectors. 
Detected bursts had shorter durations for the firing-rate-threshold and rank 
surprise detectors. This was especially noticeable for cases with fewer 
short-duration small bursts in (B) and (C) (right column and D). Arrows in (B) 
correspond to examples of a large and a small burst in (D). (D) Plotted in (D) 
are durations of identified network bursts using each method (colored bars; 
top) and identified first-stage single-channel bursts using the IS! and Rank 
surprise methods (colored raster plots; each dot is a spike). 



to assign significance values to each burst. This is indeed a desir- 
able feature if detection performs adequately. However, underesti- 
mating the number and duration of bursts appears to be inherent 
to the algorithm: The underlying Poisson or rank distribution, 
or what is "unsurprising," is created including the "surprising" 
spikes that are in bursts. Therefore, spikes in the tails of bursts, 
where ISIs are longer than in the burst core, wUl no longer be "sur- 
prising." In Figures 6, 7D, the tails of bursts are clearly excluded. 
This effect will be amplified for preparations with large percent- 
ages of spikes occurring within bursts, such as in vitro networks 



(90 ± 12% of spikes were within network bursts according to the 
ISIa;= 10 -threshold detector; mean ± SD, 9 cultures). Lowering 
the default threshold for what would be surprising entailed the 
detection of more spikes in bursts. However, the smallest bursts 
were still not detected, and, at the same time, wide stretches of 
tonic activity were now erroneously identified as bursts. Perhaps 
improved performance could be achieved by considering more 
than one distribution, with one representing a bursting regime, 
and then measure "surprise" with respect to modeled ISI distri- 
butions of the non-bursting regimes. The durations of the bursts 
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FIGURE 8 I Rate-thresholds and burst detection for data presented in 
Figure 2. The probability distributions of (A) the total number of spikes or 
(B) the total number of electrodes that detected a spike within a time 
window between 5 and 50 ms in duration (lines) are plotted. Elevated firing 



during network bursting corresponds to higher spike or electrode counts, and 
the large arrows indicate the rate-threshold for burst detection used in (C). 
(C) Detector performance for a segment of network activity (black dots). 
Colored bars indicate detected bursts for each burst detector. 



detected with the Rank Surprise and rate-threshold detectors were 
35 ± 15% and 49 ± 15%, respectively, of the durations of the 
same bursts identified by the ISI]v=io-threshold detector (mean 
± SD for 9 cultures). 

In their most basic form, rate-threshold and ISIjvf-threshold 
(including ISI-threshold) methods offer two sides of the same 
coin. The former measures the number of spikes detected within 
a fixed window of time while the latter measures how much 
time elapsed for a fixed number of spikes to occur. To empha- 
size the point, the rate-threshold detector in Figure 6 is related 
to an ISIat -threshold detector having an ISIn=8 threshold equal 
to 50 ms. The choice of number of spikes or time window as a 
threshold is arbitrary, but the former creates artificial time points 
for burst boundaries at threshold crossings, while the latter allows 



the assignment of specific time points for the first and last spike 
in putative bursts. 

The "discharge density" histograms (i.e., probability distribu- 
tions of firing rates) (Kaneoke and Vitek, 1996; Ji and Wilson, 
2007) provide a means to flip the coin and automatically deter- 
mine both thresholds. By using a fixed firing-rate histogram time 
window, the threshold N (number of spikes) was selected at the 
valley separating high and low firing rate peaks in the discharge 
density plot (c.f Figures). Flipping the coin, the method pre- 
sented here can use the N found above to select an ISI^ threshold 
(c.f Figure 2). In turn, the ISIn threshold can update the time 
window of a new discharge density histogram in order to update 
N, thus continuing an iterative process. For our data, N and 
the ISIjv threshold converged, but whether or not these values 
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are optimal for burst detection is debatable. For example for 
data presented in Figures 1 and 6, the values converged to an 
ISIj/=260 threshold equal to 200 ms and an ISIj/=i7 threshold 
equal to 220 ms, respectively. Detection appeared to work well. 
However, the small bursts in Figure 1 (i.e., less than N = 260) 
could no longer be detected, and longer ISIjv thresholds led to 
some merged bursts. For data in Figure 1, the process appeared to 
find the border between large and small bursts (c.£ compare N = 
260 to Figures). Alternatively, the simultaneous use of multiple 
pairs of N and ISIn thresholds is expected to improve detection 
performance. A simple example is given in the next paragraph, 
where an additional ISI]v=2 threshold helped to reduce the num- 
ber of merged bursts. As a last note, the valleys in the probability 
distribution histograms provided an intuitive choice for a bor- 
der but are not necessarily optimal. Finding optimal borders, 
and similarly precise measures of performance, requires having 
a ground truth about when bursts occur, which, unfortunately, 
does not exist. 

Tonic activity can merge neighboring bursts for ISI^r- 
and ISI-threshold detectors. Excluding tonic channels in a 
pre-processing stage, as discussed in the Results, offers one 
solution. Alternatively, merged bursts could be separated in a 
post-processing stage by finding time points of lowest firing 
between putative bursts as described by (Wagenaar et al, 2005; 
Pasquale et al., 2010); this technique is inherently done by rate- 
threshold detectors. However, whether or not putative bursts are 
distinct or different phases of the same burst was not always 
apparent. A further alternative that avoids pre- or post-processing 
steps is to introduce a second parameter during burst detection. 
For example, an ISIjv=2 threshold of 20 ms (50 Hz) effectively 
broke into pieces the tonic spike trains, which did not sustain 
50 Hz firing rates. Since an ISIjv=2 threshold may be difficult to 
identify (c.£ Figure 2A), such a threshold may need to be adjusted 
in an ad-hoc manner for different preparations. 

CONCLUSION 

A network burst was defined as a period of elevated firing within 
or across multiple channels. Creating a network burst detec- 
tor then consists of finding appropriate thresholds to separate 
periods of low and high firing states. We proposed an ISI-- 
threshold burst detector, where ISIn is the inter-spike-interval 
between every Nth spike in the network. The ISIn -threshold burst 
detector performed well with our data, but we acknowledge that 
other detectors may be better suited for other types of data. For 
example, the ISI]v-threshold burst detector struggled most when 
multiple channels exhibited intermittent and irregular activity. 
For single neuron bursts, an ISI-threshold burst detector (ISIjv=2) 
sufficiently distinguishes bursting regimes (Pasquale et al., 2010). 

Our approach offered a number of useful features including: 
a simple and computationally efficient implementation, no need 
for ad-hoc or post-hoc criteria, and precise assignment of burst 
boundary time points. The choice of N is not obvious, however, 
and could be considered to be an ad-hoc decision based on famil- 
iarity with the structure of the spike trains. We chose N = 10 as 
it was large enough to demonstrate a deepening of valleys in the 
probability distributions (Figure 2A) and small enough to allow 
a wide range of burst sizes to be detected. Choosing N equal to 



the minimum burst size would give the least sensitivity to tonic 
or noisy channels, but a minimum burst size may not be obvi- 
ous (c.f Figures). On the other hand, choosing the smallest N 
that provides easily separable peaks in the ISIjv probability distri- 
bution would not bias detection toward larger bursts, allowing 
identification and analysis of a greater range of neuronal and 
network dynamics. 
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SUPPLEMENTARY MATERIAL 

The Supplementary Material for this article can be found 
online at: http://www.frontiersin.org/Computational_ 
Neuroscience/10.3389/fncom.2013.00193/abstract 

Supplemental Figure 1 | Rate-thresholds and burst detection for 
spike-sorted data presented in Figures 2 and 8. The probability 
distributions of (A) the total number of neurons spiking or (B) the total 
number of electrodes that detected a spike within a time window 
between 5 and 50 ms in duration (lines) are plotted. Elevated firing during 
network bursting corresponds to higher neuron or electrode counts, and 
the large arrows indicate the rate-threshold for burst detection used in (C) 
and (D). Hundred and two electrodes detected spikes, and from these, 62 
individual neurons were manually identified (spike-sorted) based on having 
distinct spatio-temporal activity profiles (Franke et al., 2012a). (C) 
Histograms of neuron or electrode counts for 5 or 50 ms time windows 
for the network activity presented in (D). The rate-thresholds detected in 
(A) and (B) are plotted as dotted lines, and a burst is detected whenever a 
count exceeds the rate-threshold. (D) Detector performance for a 
segment of network activity (black dots). Colored bars indicate detected 
bursts for each burst detector. 

Supplementary Code 1 1 Matlab code for creating ISI/v histogram plots in 
order to choose an ISI/v threshold. 

Supplementary Code II | Matlab code for ISI/v burst detection. 
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